
# library(pacman)
# p_load(dplyr, readr, fixest, modelsummary, kableExtra) 
library(dplyr)
library(readr) 
library(fixest)
library(modelsummary)
library(kableExtra) 

options(modelsummary_factory_default = 'kableExtra')

alldata1 <- read_rds("replication_data1.rds")
alldata2 <- read_rds("replication_data2.rds")
alldata3 <- read_rds("replication_data3.rds")

# In the paper, only the main independent variable is shown. 
cm <- c("weight" = "ln(TC connections)") 

mod_all0 <- feols(logtrade ~ weight | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all1 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year,
                  cluster = "dyad_id+year", data = alldata1)

mod_all2 <- feols(logtrade ~ weight + patents_cty_1 + patents_cty_2 +
                    fta_wto + wto_dyad | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata1)

mod_all3 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur | 
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata2)

mod_all4 <- feols(logtrade ~ weight + 
                    fta_wto + wto_dyad +
                    dem_dyad + pta + comcur +
                    alliance + rivalry |
                    dyad_id+cty_1+cty_2+year, 
                  cluster = "dyad_id+year", data = alldata3)

models <- list(mod_all0, mod_all1, mod_all2, mod_all3, mod_all4)
names(models) <- c("Baseline", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++")

rows <- tibble::tribble(~start, ~mod_all0, ~mod_all1, ~mod_all2, ~mod_all3, ~mod_all4, 
                        "Controls", "No", "Gravity", "Gravity+R&D", "Gravity+", "Gravity++",
                        "Time series", "2004-2022", "2004-2021",  "2004-2021", "2004-2015", "2004-2011") 

tabsummary <- modelsummary(models, 
                           fmt = 3,
                           coef_map = cm,
                           #output = "latex", # This line produces latex output of the table
                           stars = TRUE,
                           gof_omit = 'AIC|BIC|Within|Std.Errors|R2|FE',
                           notes = c("Fixed effects by dyad, country and year, clustered standard errors by dyad and year.",
                                     "Zero imputation on dyads with missing on TC connections.",
                                     "Gravity controls: GDP, population, distance between capitals, common language, regional trade agreement, WTO dyad.",
                                     "Gravity+R&D controls: Adds to Gravity patents per country.",
                                     "Gravity+ controls: Adds to Gravity+ democratic dyad, preferential trade agreement, common currency.",
                                     "Gravity++ controls: Adds to Gravity++ strategic rivalry, alliance."),
                           add_rows = rows)

tabsummary %>%
  add_header_above(c(" " = 1, "Dependent variable: ln(Dyadic trade) (UN Comtrade)" = 5)) %>%
  kable_styling(font_size = 8, full_width = F) 
